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A theoretical model is presented to compute the electronic band structures and optical absorption 
spectra for twisted incommensurate few-layers graphene (tLFG) systems of arbitrary architecture. 

This is accomplished using an integrated virtual crystal approximation and tight-binding (VCA- 
TB) method, where the VGA is achieved by a mathematical averaging formalism developed over 
the quasi-infinite ensemble of interlayer bond configurations. The results show that the low-energy 
electronic band structures of the incommensurate tLFG systems of N graphene layers are formed of N 
overlapping Dirac cones centered on the K-points of the Brillouin zones of the graphene layers. This 
yields effective gaps between the saddle points of the valence and conduction bands of these tFLG 
systems, that are tunable with the incommensurate twist angles. The optical absorption spectra are 
also calculated for these systems; in particular, they display clear spectral peaks which correspond 
to the accessible electronic transitions across the gaps between the valence and conduction bands. 

The optical absorption results highlight the potential of these incommensurate tFLG systems for 
applications in graphene based tandem cells. 
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I. INTRODUCTION 

Bilayer graphene may be stacked in three different con¬ 
figurations, namely the AB-type (Bernal stacking), AA- 
type and twisted-type. The twisted bilayer graphene 
(tBLG) is the most general, consisting of two sheets of 
graphene rotated by an arbitrary twist angle 9 other than 
an integer multiples of 60°. Each of the three stacking 
types exhibits a different low energy electronic structure, 
and hence a different optical response mm- The inter¬ 
layer interaction in tBLG, in particular, modifies signif¬ 
icantly the low-energy band structure and yields novel 
electronic features not present in monolayer graphene or 
other bilayer graphenes [1]-[I3]- The tBLG, therefore, 
has the richest optical absorption spectrum, owing to its 
complex low energy band structure, where the optical re¬ 
sponse is observed to be highly tunable through the twist 
angle 9 [I1].[T5]. 

Gommensuration in tBLG occurs at very specific twist 
angles m- As a consequence, experimentally pre¬ 
pared tBLG structures are in general incommensurate. 
The electronic band structures of bilayer and few-layers 
graphene have been investigated intensively using ab¬ 
sorption measurement experiments m-m- Theoreti¬ 
cally, various models have been developed to calculate 
the electronic structure in commensurate tBLG 
[IS],|1S]-I53]. These models were also used to calculate 
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the optical absorption for tBLG for a wide range of ro¬ 
tation angles, doping conditions, and interlayer poten¬ 
tial asymmetry MM, Ea-EH]. However, the optical 
properties of incommensurate twisted bilayer (tBLG) and 
few-layers graphene (tFLG) have not been investigated 
yet on the basis of a systematic theory for incommensu¬ 
rate heterostructures. 

Incommensurate tBLG constitutes a disordered het¬ 
erostructure characterized by a quasi-continuum set of 
interlayer bonds, and a quasi-infinite set of configura¬ 
tions for carbon atoms at the bilayer interface. In a pre¬ 
vious work m, we have developed a non-trivial Virtual 
Grystal Approximation - Tight Binding (VCA-TB) the¬ 
ory to study the electronic properties in such systems. 
Our VGA-TB theoretical approach determines the VGA 
effective medium for the incommensurate tBLG system, 
derived in a direct manner from the quasi-infinite set of 
structural configurations at the bilayer interface. The 
theory consequently determines the VGA-TB Hamilto¬ 
nian for the effective VGA system and yields the tBLGs 
electronic structure and density of states for any incom¬ 
mensurate twist angle 9. 

The purpose of the present work is to extend the VCA- 
TB theoretical approach to study the electronic structure 
and the optical absorption properties in incommensurate 
twisted few-layers graphene (tFLG) with arbitrary stack¬ 
ing geometries. The band structure for incommensurate 
tFLG calculated by the VCA-TB method is found to be 
formed of overlapped Dirac cones, equal in number to 
that of the graphene sheets forming the incommensurate 
heterostructure. The overlap between these Dirac cones 
gives rise to numerous saddle points in the low-energy 
band structure and consequently to van Hove singular¬ 
ities (VHS) in the density of states (DOS). The optical 
absorption spectrum, calculated in the framework of the 
VCA-TB theory, is found to be characterized by a series 
of peaks associated with the angle dependent van Hove 
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singularities. The transition energies corresponding to 
optical absorption peaks are hence directly tunable via 
the rotation angles defining the twists in the stacking 
configurations of the incommensurate tFLG. 

The rest of this paper is organized as follows. In Sec¬ 
tion II we extend the VCA-TB model and apply it to 
determine the band structures and densities of states for 
incommensurate twisted tFLG. In Section III, we calcu¬ 
late the optical absorption for tFLG for various incom¬ 
mensurate architectures and analyze the results in view 
of the characteristic features of the band structure. The 
conclusions are discussed in Section IV. 
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II. BAND STRUCTURE FOR 
INCOMMENSURATE TWISTED FEW-LAYERS 
GRAPHENE 

Consider an incommensurate heterostructure of N 
twisted graphene layers stacked on top of each other. The 
layers are labeled by n=l, 2,.., N. Any layer n -I- I is as¬ 
sumed twisted with respect to adjacent layer n by an 
incommensurate twist angle 9n,n+i- The A and B sub¬ 
lattices in the graphene layer n are denoted by A„ and 
Bn, and the corresponding atoms are labeled as A„(s) 
and Bn{s)\ Fig.I shows only the first eighteen sites of 
the infinite set of sites per layer. 

In the VCA-TB model, we consider a referential 
couple of sites in each graphene layer n, denoted by 
A„(0), Bn{0). The incommensurate tFLG is hence repre¬ 
sented by a basis of 2N irreducible atomic sites, namely 
Ai(0),Bi(0),.., Ajv(0),B7v( 0). In the VCA-TB method, 
the effective VGA medium is derived as a direct system 
average from the quasi-infinite set of structural configura¬ 
tions at bilayer interfaces within the tFLG. The effective 
medium is hence periodic and one may write the total 
electronic wave function for the tFLG as 

N N 

n—1 n—1 

In our formalism, we consider a single p^-orbital per site. 
The wave function 'I'k;a(r) = l/^/W — 

fa), a = An or Bn, is the VCA-TB Bloch wave function 
from the a sublattice in layer n. <I>q, denotes a 2 p 2 atomic 
orbital and the letter V is the number of unit cells in 
the crystal, is the position of an a type atom with 
respect to a chosen origin. The coefficient Ca (a = A„, 
or Bn) is the a sublattice contribution to the total wave 
function and k is the wave vector. The band structure 
in the incommensurate tFLG system may be determined 
by solving the secular equation 
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FIG. 1: Figs, a) and b) present schematic representations of 
the An and Bn snblattices in layer n of the incommensurate 
tFLG systems, with representative sublattice atomic sites. 


///d3r4'*,^(r)5'k;/3(r) and = (|^'k;a|-ff|^'k;/3) = 

///d^i’^k;aW-ff^k;/3(r); a, ^ G = I,...,iV}. 

The in-plane entities of the VCA-TB Hamiltonian ma¬ 
trix H may be readily written as |39] 


HAnAn — HSnBn “ 


2p 


I + 2si ^ cos(ai(n).k) 


i=l 


27iy^cos(a^(7t).k), 


HAnBn — 

= £’ 2 p |so7k(n)-b S 2 2cos(ai(n).k)-b I 

+ 7o7k(7i)-b72 2cos(ai(n).k)-b . 


det[H - E{k)S]=0, (I) 

or equivalently by determining the eigen-energies A(k) 
of the matrix S~^H. The matrices H and S respec¬ 
tively denote the VCA-TB Hamiltonian and overlap 
matrices with matrix entities Sap = ('I'k;a|'I'k;/ 3 ) = 


Here we have defined the lattice vectors SLi{n) for 
j=I, 2, 3, 4 for the graphene layer n as aiin) = 
R (^* 1,2 + ^ 2,3 + + dn-i,n) ^^(1), where R denotes the 2 x 
2 rotation matrix, and the vectors ai(l) = a(-\/3,0), 
a 2 (l) = ^a(l,V3), a 3 (I) = a 2 (l) - ai(l), a 4 (I) = 
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a(0, 2). The parameters 70 , 71,72 and so,si,S 2 denote 
the in-plane first, second and third nearest neighbor hop¬ 
ping and overlap integrals, respectively. E 2 p denotes the 
onsite energy and 7k(?^) = 1 + -I- e*‘^ 2 (n) k^ 

Due to the incommensurate nature of the tFLG stack¬ 
ing, the carbon atoms in these systems are characterized 
by a quasi-infinite ensemble of structural configurations 
with a quasi-continuum set of interlayer bond lengths and 
orientations. The key idea in the VCA-TB approach is 
to determine the ensemble average for the interlayer hop¬ 
ping integrals between neigboring sites, over the quasi¬ 
infinite ensemble of structural configurations. We use 
the environment-dependent tight binding approach, m, 
[40j . to compute the values of the interlayer hopping in¬ 
tegrals and of their upper limit as a function of the ran¬ 
domly distributed bond lengths and configurations. The 
VGA is achieved by a mathematical averaging formalism 
developed over the quasi-infinite ensemble of interlayer 
bond configurations |41j . This model approach has been 
applied successfully for the tBLG systems with arbitrary 
incommensurate twist angles m- We compute in the 
present work for the generalized tFLG systems, the en¬ 
semble averages of the interlayer elements of the VGA-TB 
Hamiltonian matrix H as 

Ha^a^ = [t(A„(0)A™(0)) + t(A„(0)A™(l)) 
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X ^cos((r^^(s)-rA^(o)) -k)] 

S^l 

^ (^n,m —1 T • 

ryi„(s) denotes the position vector of the A„(s) atom as in 
Fig.l. The quantity t (A„(0)Am(s)) denotes the average 
for the interlayer hopping integral between atoms A„(0) 
and Ajn{s), computed using the VGA-TB approach m- 
The Kronecker deltas represent the interactions between 
adjacent graphene layers. 

Next we calculate the elements of the overlap S ma¬ 
trix. The interlayer overlap integrals are neglected owing 
to the large perpendicular distance d± between the two 
graphene layers [dj_ « 3.35A). Note in this case that 
the S matrix eelements in the plane of a graphene layer 
are invariant under the twist angle. The nonzero matrix 
elements of the overlap matrix S may be represented as 

3 

= Sb„b„ = 1 + 2si ^ cos(aAn).k), 

2 = 1 


= Sb^a^ = ■So7k(ri)+S2 [2cos(ai(n).k) -b 


At this stage, one can solve numerically the secular 
Eq.(l) and determine the band structure for any given 
tFLG system with any number of graphene sheets iV, 
and incommensurate twist angles 9n,n+i- The numerical 
values for the various averaged interlayer hopping 
integrals are presented in reference |13] . whereas the 
numerical values for the onsite energy A 2 p, the in-plane 
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FIG. 2: The numerically calculated low-energy electronic 
band structure for the incommensurate tFLG {N = 4,0i2 = 
2 °,023 = 2°,034 = 2°), plotted over the juxtaposed set of 
the Brillouin zones corresponding to the incommensurably 
twisted graphene layers forming the graphene heterostructure. 
Our numerical results show that the band structure of a tFLG 
composed of N graphene layerss is formed of N overlapped 
Dirac cones. 



hopping integrals 70 , 71 , 72 , and the overlap integrals 
so,si,S 2 are taken from reference |39) . 

To illustrate the numerical results, we present in Fig.2 
the calculated low energy electronic band structure for 
the tFLG with four incommensurable graphene sheets 
(^N = 4,012 = 2 °,023 = 2°,034 = 2°); the results demon¬ 
strate that the band structure of a tFLG composed of N 
graphene sheets is formed of N overlapped Dirac cones. 
This overlap significantly gives rise to numerous saddle 
points for the valence and conduction bands in the low 
energy electronic band structure for these systems. 

To further illustrate the behavior of the valence and 
conduction bands, and their saddle points, with the 
variation of the graphene incommensurate stacking, we 
present in Figs.3 the low energy electronic band struc¬ 
tures for two tBLG and three tFLG systems, along the 
high symmetry axes TKi and KiMi in the Brillouin 
zone of layer 1. To zoom on the band structures near the 
saddle points, we opt in the BZ without loss of generality 
for the points Gi and A^i along TKi and KiMi respec¬ 
tively, with FGi = O.SSFATi and KiNi = 0.33KiMi. 

In particular. Figs.3a, 3b, and 3b, correspond to the 
band structures for two tBLG systems with increasing 
twist angles, namely {N = 2, 612 = 1°) and {N = 2,0i2 = 
2°), whereas Figs.3c, 3d and 3e, correspond to the band 
structures for three tFLG systems with increasing twist 
angles, namely {N = 3,0i2 = 1°,023 = 1°), (^ = 3,0i2 = 
2°,023 = 2°), and (N = AA2 = 1“,023 = 1“,034 = 1°), 
respectively. Note that the Fermi energy is taken as a 
base line reference at Ep = 0. 

The saddle points in the low energy electronic band 
structures of Figs.3 yield the van Hove singularities in 
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FIG. 3: Calculated low energy electronic band structures for 
tBLG and tFLG systems, along the high symmetry axes FKi 
and K\M\. Plotted over the BZ of a reference graphene 
layer in a given system, these correspond to: (a) tBLG 
{N = 2fiu = 1°), (b) tBLG {N = 2,6^ = 2°), (c) tFLG 
{N = 3,6)i 2 = 1°,023 = 1°), (d) tFLG (TV = 3,0i2 = 2°,023 = 
2°), and (e) tFLG [N = 4,0i2 = 1°,023 = 1°,034 = 1°). The 
results show the saddle points of the valence and conduction 
bands with associated gaps due to the overlapping of the Dirac 
cones from the system graphene layers. 


the corresponding densities of states (DOS), calculated 
along TKi and KiMi and presented in Figs.4. It is 
clear that the positions of the spectral singularities may 
be directly tuned via the incommensurate twist angles 
of the tFLG systems. Note that the arrows in Figs.4 
highlight the electronic transitions by optical absorption 
across the gaps between the saddle points of the valence 
and the conduction bands. 


III. OPTICAL ABSORPTION 


€ is a positive infinitesimal and lo = E/Ti. The con¬ 
stants 5s and gv denote the spin and electronic valley 
degeneracies, respectively {gs = g^ = 2). The opera¬ 
tor jx is the current operator defined in our formalism as 
jx = —dS~^H/dkx. The function / is the Fermi function 
and la, k > is the a sub-band eigenstate of the operator 

S~^H with momentum k and eigen-energy i?a(k). 

Given the calculated electronic states within the 
VGA-TB approach, we use Eq.(2) to calculate the 
dynamic conductivity <Jxx{E) for the tFLG systems 
along the high symmetry directions TKi and KiMi. 
This calculation yields the optical absorption spectra 
Re(Jxx{E)/Reamono{E) for these tFLG systems, normal¬ 
ized with respect to the optical absorption for a graphene 
monolayer along the same symmetry directions. 

In Figs.5, we plot the normalized optical absorption 
spectra at room temperature over the low energy inter¬ 
vals of the band structures in the neighborhood of the 
Fermi level Ep = 0. The series of results of Figs.5 corre¬ 
spond to the series of the systems of Figs. 3 and 4. 

There are a few general remarks which can be made 
about Figs.5. The normalized optical absorption for the 
tFLG systems in the neighborhood of the Fermi energy 
is observed to approach unity. This is a direct conse¬ 
quence of the fact that the energy spectrum of the tFLG 
systems becomes identical to that of monolayer graphene 
very close to the Ki Dirac point as in Figs. 3. The cal¬ 
culated spectra further display peaks which correspond 
to the electronic transitions from the saddle points in 
the valence bands to the saddle points in the conduction 
bands. The peaks in the optical absorption spectra are 
highlighted using arrows as in Figs.4, in order to identify 
each peak with the corresponding transition. Our numer¬ 
ical results demonstrate the possibility to tune the tran¬ 
sition energy corresponding to these peaks by the manip¬ 
ulation of the incommensurate twist angles of the layered 
heterostructures. Finally, the calculated normalized ab¬ 
sorption spectra at relatively high energies are observed 
to converge to a value equal to the number of graphene 
layers N. This result follows from the fact that the elec¬ 
tronic structures at such energies are similar to those for 
monolayer graphene but with a number of modes N times 
greater than in monolayer graphene. 


The optical absorption of light incident along the nor¬ 
mal to the layers of the tFLG system, is defined as the 
real part of the dynamic conductivity a^x, where the x- 
axis is in the plane of the layers. The dynamic conductiv¬ 
ity may be calculated then along any given direction A 
in the reciprocal space of the BZ for any of the graphene 
sheet layers, via the Kubo formula 

.pN _ . 9s9ve'^ f f{Ea{k)-f{Ei3{k)) 

^ ^uj + ieJ^27T^^E-Epik) + E^(k)+ie 

X (a;,k|ja,|/3,k^ {^^,k|j2;|a,k^ . (2) 


IV. CONCLUSIONS 

We have developed a theoretical approach to model 
and compute the electronic band structures and opti¬ 
cal absorption spectra for twisted incommensurate few- 
layers graphene tLFG systems of arbitrary architec¬ 
ture. This is accomplished using the integrated VGA- 
TB method and the environment-dependent tight bind¬ 
ing approach, where the VGA is achieved by a mathemat¬ 
ical averaging formalism developed over the quasi-infinite 
ensemble of bond configurations. 

The numerical calculations demonstrate that the elec- 
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a) N=2, 012=1° b) N=2, 012=2° c) N=3, 0i2=l°, 023=1° 





FIG. 4: Calculated electronic DOS for the tBLG and tFLG systems, detailed in Figs.3 as (a), (b), (c), (d), and (e), along the 
high symmetry axes FTFi and KiMi\ the results are plotted over the first BZ of a reference graphene layer. 





d) N=3,012=2°, 023 =2° 



e) N=4, 012 = 1°, 023 = 1°, 034 = 1° 



FIG. 5: Calculated room temperature optical absorption spectra for the tBLG and tFLG systems, (a) and (b), and (c), (d), 
and (e) respectively, detailed in Figs.3, normalized with respect to the optical absorption spectrum for monolayer graphene. 
The absorption spectra are calculated from contributions along the high symmetry axes TKi and KiMi. 


tronic band structures of the incommensurate twisted few 
layers graphene of N sheets, are formed of N overlapping 
Dirac cones centered on the K-points in the Brillouin 
zones of the graphene sheets, generating twist-angle de¬ 
pendent saddle points in the band structures, and conse¬ 
quently effective gaps between the valence and conduc¬ 
tion bands with corresponding DOS van Hove singular¬ 


ities. The virtual crystal approximation-tight binding 
method also permits the determination of the optical ab¬ 
sorption spectra for these tFLG systems. The calculated 
optical absorption spectra are found to be highly tun¬ 
able via the twist angle. In particular, these spectra 
display clear peaks which correspond to the electronic 
transitions from the saddle points in the valence band 
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to those in the conduction band. The optical absorp¬ 
tion results highlight the potential importance of the in¬ 


commensurate tFLG systems for possible technological 
applications in graphene based tandem cells. 
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